﻿using System;
using System.Collections.Generic;
using System.Collections.Specialized;
using System.Text;
using System.Windows.Forms;
using System.Drawing;
using System.IO;

namespace GeoFly
{
    /// <summary>
    /// DEM处理，提取DEM信息和自然水系
    /// </summary>
    public class DEMRiverNet
    {
        /*======================================================================+
        +							DEM提取水系类定义开始						+
        +=======================================================================+
        +		 class name:	DEMRiverNet										+
        +		 Function:		(1)填洼处理										+
        +						(2)坡度、坡向计算	(2005-03-28:通过检验)		+
        +						(3)水流流向计算			(2005-03-28)			+	
        +						(4)水流累积矩阵		(2005-03-29:通过检验)		+
        +						(5)流域河网提取		(2005-03-30:通过检验)		+
        +						(6)流域边界提取		(2005-03-30:通过检验)		+
        +						(7)												+
        +						(8)												+
        +=======================================================================*/
        public DEMRiverNet(GridLayer DemData)
        {
            //this.m_nodata = -9999;
            this.demgrid = DemData;            
        }

        
        /*======================================================+
        +			   地形、河道参数计算函数定义				+
        +======================================================*/

        /// <summary>
        /// 用Horn方法计算坡度
        /// </summary>
        /// <returns></returns>
        public double[,] CalcSlope()
        {
            //生成存放坡度的矩阵
            double[,] result = new double[demgrid.rowCount, demgrid.colCount];
            double[,] v = demgrid.Values;//用v来代替demgrid.pGridValue变量

            for (int i = 1; i < demgrid.rowCount - 1; i++)
            {
                for (int j = 1; j < demgrid.colCount - 1; j++)
                {
                    if (v[i, j] != demgrid.nodata)
                    {
                        //开始计算坡度
                        double p = (v[i - 1, j - 1] + 2 * v[i - 1, j] + v[i - 1, j + 1]) - (v[i + 1, j - 1] + 2 * v[i + 1, j] + v[i + 1, j + 1]);
                        double q = (v[i - 1, j - 1] + 2 * v[i, j - 1] + v[i + 1, j - 1]) - (v[i - 1, j + 1] + 2 * v[i, j + 1] + v[i + 1, j + 1]);
                        result[i, j] = Math.Sqrt(p * p + q * q) / 8 / demgrid.resolution;
                        //注意上面用的公式与张东的不同
                    }
                }
            }
            this.m_pSlope = result;
            return result;
        }
        /// <summary>
        /// 用Horn方法计算坡向
        /// </summary>
        /// <returns></returns>
        public double[,] CalcAspect()
        {
            //生成存放坡向的矩阵
            double[,] result = new double[demgrid.rowCount, demgrid.colCount];
            double[,] v = demgrid.Values;//用v来代替demgrid.pGridValue变量
            for (int i = 1; i < demgrid.rowCount - 1; i++)
            {
                for (int j = 1; j < demgrid.colCount - 1; j++)
                {
                    if (v[i, j] != demgrid.nodata)
                    {                        
                        double p = (v[i - 1, j - 1] + 2 * v[i - 1, j] + v[i - 1, j + 1]) - (v[i + 1, j - 1] + 2 * v[i + 1, j] + v[i + 1, j + 1]);
                        double q = (v[i - 1, j - 1] + 2 * v[i, j - 1] + v[i + 1, j - 1]) - (v[i - 1, j + 1] + 2 * v[i, j + 1] + v[i + 1, j + 1]);
                        result[i, j] = Math.Atan(q / p);
                    }
                }
            }
            this.m_pAspect = result;
            return result;
        }
        /// <summary>
        /// 在DEMlayer矩阵中找出属于水体的像元，连片的像元认为构成了一个水体
        /// </summary>
        /// <param name="DEMlayer"></param>
        /// <returns></returns>
        public List<WaterBody> FindWaterBody(GridLayer DEMlayer, int[,] DirectionMatrix)
        {
            List<WaterBody> waterBodies = new List<WaterBody>();
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在搜索水体";
            bar.Show();
            int rowCount = this.demgrid.rowCount;
            int colCount = this.demgrid.colCount;
            for (int rowIndex = 0; rowIndex < rowCount; rowIndex++)
            {
                for (int colIndex = 0; colIndex < colCount; colIndex++)
                {
                    if (DirectionMatrix[rowIndex, colIndex] >0)
                        continue;
                    Cell tempCell=new Cell(rowIndex,colIndex);
                    //如果在这些水体中的任何一个中发现了本像元，就跳过
                    bool finded = false;
                    foreach (WaterBody body in waterBodies)
                    {
                        if (body.WaterCells.ContainsKey(tempCell))
                        {
                            finded = true;
                            break;
                        }
                    }
                    //如果该像元未进行任何水体，则搜索它是否属于一个水体的一部分
                    if (!finded)
                    {
                        WaterBody mbody = this.SearchWaterBody(DEMlayer, tempCell);
                        if (mbody != null)
                        {
                            waterBodies.Add(mbody);
                        }
                    }
                }
                bar.progressBar1.Value = (int)(rowIndex * 100.0 / DEMlayer.rowCount);
            }
            bar.Close();
            return waterBodies;
        }
        /// <summary>
        /// 从像元cell开始，以漫水法搜索高程值相同的邻接像元，把这些高程值相同的连片像元组成水体。
        /// </summary>
        /// <param name="DEMlayer"></param>
        /// <param name="cell"></param>
        /// <param name="waterBodies"></param>
        /// <returns></returns>
        private WaterBody SearchWaterBody(GridLayer DEMlayer, Cell cell)
        {
            double[,] Values = DEMlayer.Values;
            double ReferValue = Values[cell.row, cell.col];
            Stack<Cell> stack = new Stack<Cell>();
            WaterBody body = new WaterBody();
            stack.Push(cell);
            body.WaterCells.Add(cell,-1);
            while (stack.Count > 0)
            {
                Cell tempCell = stack.Pop();
                
                int i = tempCell.row;
                int j = tempCell.col;
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    int I = i + index / 3 - 1;
                    int J = j + index % 3 - 1;
                    if (I < 0 || I >= DEMlayer.rowCount)
                        continue;
                    if (J < 0 || J >= DEMlayer.colCount)
                        continue;
                    double value=DEMlayer.Values[I, J];
                    Cell newCell = new Cell(I, J);
                    if (!body.WaterCells.ContainsKey(newCell) && value == ReferValue)
                    {
                        stack.Push(newCell);
                        body.WaterCells.Add(newCell, -1);
                    }
                }
                //if (stack.Count > 1000)
                //    throw new Exception("");
            }
            if (body.WaterCells.Count > 1)
                return body;
            else
                return null;

        }
        private int[,] marks;
        private int ProcessCount=0;
        public int[,] marks2;
        
        /// <summary>
        /// 返回marks的值
        /// </summary>
        /// <param name="row"></param>
        /// <param name="col"></param>
        /// <returns></returns>
        public int Marks(int row,int col)
        {
            if (row < 0 || row>= this.demgrid.rowCount || col< 0 || col>= this.demgrid.colCount)
                return 0;
            else
                return marks[row, col];
        }
        private void Enqueue(SortedDictionary<double, Queue<Cell>> GQueue, Cell cell, double height)
        {
        	 if (GQueue.ContainsKey(height))
                {
                    GQueue[height].Enqueue(cell);
                }
                else
                {
                    Queue<Cell> queue = new Queue<Cell>();
                    queue.Enqueue(cell);
                    GQueue[height]=queue;
                }
        }
        
        /// <summary>
        /// 填洼算法
        /// 设定的方向编码方案如下：
        ///     i -------------
        ///     0 | 0 | 1 | 2 |
        ///       -------------
        ///     1 | 3 | 4 | 5 |
        ///      -------------
        ///     2 | 6 | 7 | 8 |
        ///       -------------
        ///         0   1   2 j
        ///  若当前像素(i,j)的流向编码为m，则此编码指向的相邻相素的行标和
        ///  列标分别为
        ///         I=i+m/3-1
        ///         J=j+m%3-1
        ///  当前像素某方向上的相邻像素的编码为a，则本像素位于相邻像素的方
        ///  向编码为8-a。其中编码4表示中心位置，在流向编码中应该避免使用
        ///  该码
        /// </summary>
        /// <param name="DirectionMatrix">方向矩阵</param>
        /// <param name="dem">数字高程模型</param>
        /// <returns></returns>
        
        public void CheckPit_LiuYonghe2()
        {
        	int DequeCount=0;
            Random random = new Random();
            //创建一个与流域矩阵大小相同的矩阵，用于标记流域是否是水体，是水体为0，不是水体为1；默认不是水体
            marks = new int[demgrid.rowCount, demgrid.colCount];
            for (int row = 0; row < demgrid.rowCount; row++)
            {
                for (int col = 0; col < demgrid.colCount; col++)
                {
                    marks[row, col] = 1;
                }
            }
          
            marks2=new int[demgrid.rowCount,demgrid.colCount];
            //创建优先队列对象
            SortedDictionary<double, Queue<Cell>> GQueue = new SortedDictionary<double, Queue<Cell>>();
            for (int rowIndex = 0; rowIndex < demgrid.rowCount; rowIndex++)
            {
                //第一列的元素
                double tempValue1 = demgrid.Values[rowIndex, 0];
                //最后一列的元素
                double tempValue2 = demgrid.Values[rowIndex, demgrid.colCount - 1];
                Cell cell1 = new Cell(rowIndex, 0);
                this.Enqueue(GQueue,cell1,tempValue1);

                Cell cell2 = new Cell(rowIndex, demgrid.colCount - 1);
                this.Enqueue(GQueue,cell2,tempValue2);
                
            }
            for (int colIndex = 0; colIndex < demgrid.colCount; colIndex++)
            {
                double tempValue1 = this.demgrid.Values[0, colIndex];
                double tempValue2 = this.demgrid.Values[demgrid.rowCount - 1, colIndex];
                Cell cell1 = new Cell(0, colIndex);
                this.Enqueue(GQueue,cell1,tempValue1);
                Cell cell2 = new Cell(demgrid.rowCount - 1, colIndex);
                this.Enqueue(GQueue,cell2,tempValue2);
            }

            //ProgressBar bar = new ProgressBar();
            //bar.Text = "正在填洼";
            //bar.Show();
            //int count = 0;

            while (GQueue.Count > 0)
            {
                SortedDictionary<double, Queue<Cell>>.Enumerator en = GQueue.GetEnumerator();
                en.MoveNext();
                double minValue = en.Current.Key;
                Queue<Cell> queue = GQueue[minValue];//获取最小高程的一组像素的队列
                
                double increment = 0;
                while (queue.Count > 0)
                {
                    Cell tCell = queue.Dequeue();
                    DequeCount+=1;
                    int i = tCell.row;
                    int j = tCell.col;
                    marks2[i,j]=DequeCount;
                    //if (demgrid.Values[i, j] == 3)
                    //{
                    //}
                    if (marks[i, j] > 0)
                        marks[i, j] = 0;
                    else
                        continue;
                    //count++;
                    //if (count % 100 == 0)
                    //    bar.progressBar1.Value = (int)(count * 1.0 / demgrid.rowCount / demgrid.colCount * 100);
                    List<int> set = new List<int>();
                    for (int n = 0; n < 9; n++)
                    {
                        set.Add(n);
                    }
                    for (int n = 0; n < 9; n++)
                    {
                        int rnd = random.Next(set.Count);
                        int index = set[rnd];
                        set.RemoveAt(rnd);
                    //for(int index=0;index<9;index++)
                    //{
                   
                        if (index == 4)
                            continue;
                        int I = i + index / 3 - 1;
                        int J = j + index % 3 - 1;
                        if (I < 0 || I >= this.demgrid.rowCount)
                            continue;
                        if (J < 0 || J >= this.demgrid.colCount)
                            continue;
                        if (marks[I, J] == 0)
                            continue;
                       
                        if (demgrid.Values[I, J] <= minValue)
                        {
                            if (marks[I, J] != 2)
                            {
                                queue.Enqueue(new Cell(I, J));
                                marks[I, J] = 2;//访问但没填充过的标记
                            }
                            increment += 1e-8;
//                            if (demgrid.Values[I, J] > 0) //对于不是流域以外的海洋平面，不需要填加增值
//if(increment>0.0001)
//	throw new Exception("");
							demgrid.Values[I, J] = minValue;
//                            else
//                                demgrid.Values[I, J] = minValue;
                        }
                        else
                        {
                        	double height=demgrid.Values[I, J];
                        	this.Enqueue(GQueue,new Cell(I, J),height);
                        }
                    }
                }
                GQueue.Remove(minValue);
            }
            //bar.Close();
        }
        /// <summary>
        /// 按照planchon原始算法写的算法
        /// </summary>
        /// <param name="DirectionMatrix"></param>
        /// <param name="dem"></param>
        public void CheckPit_Planchon(int[,] DirectionMatrix, GridLayer dem)
        {
            //while (true)
            //{
                this.Init_W();
                for (int row = 0; row < this.demgrid.rowCount; row++)
                {
                    this.Dry_upward_cell(row, 0);
                    this.Dry_upward_cell(row, demgrid.colCount - 1);
                }
                for (int col = 1; col < this.demgrid.colCount - 1; col++)
                {
                    this.Dry_upward_cell(0, col);
                    this.Dry_upward_cell(demgrid.rowCount - 1, col);
                }

                for (int scans = 0; scans < 8; scans++)
                {
                    int R = R0[scans];
                    int C = C0[scans];
                    something_done = false;
                    while (true)
                    {
                        if (W[R, C] > demgrid.Values[R, C])
                        {
                            for (int index = 0; index < 9; index++)
                            {
                                if (index == 4)
                                    continue;
                                double epsilon = 1e-7;
                                if (index % 2 == 0)
                                    epsilon *= 1.4142135;
                                int I = R + index / 3 - 1;
                                int J = C + index % 3 - 1;
                                if (I < 0 || I >= this.demgrid.rowCount)
                                    continue;
                                if (J < 0 || J >= this.demgrid.colCount)
                                    continue;
                                if (demgrid.Values[R, C] > this.W[I, J] + epsilon)
                                {
                                    this.W[R, C] = demgrid.Values[R, C];
                                    something_done = true;
                                    this.Dry_upward_cell(R, C);
                                    break;
                                }
                                if (this.W[R, C] > this.W[I, J] + epsilon)
                                {
                                    this.W[R, C] = this.W[I, J] + epsilon;
                                    something_done = true;
                                }
                            }
                        }
                        Cell cell = this.NextCell(R, C, scans);
                        if (cell.row==-1)
                            break;
                        R = cell.row;
                        C = cell.col;
                    } //end while
                    if (something_done == false)
                        return;

                }//end for
                this.demgrid.Values = this.W;            
        }
        public void Dry_upward_cell(int i, int j)
        {
            depth += 1;
            if (depth > 20000)
            {
                depth--;
                return;
            }

            for (int index = 0; index < 9; index++)
            {
                if (index == 4)
                    continue;
                int I = i + index / 3 - 1;
                int J = j + index % 3 - 1;
                if (I < 0 || I >= this.demgrid.rowCount)
                    continue;
                if (J < 0 || J >= this.demgrid.colCount)
                    continue;
                double epsilon = 1e-10;
                if (index % 2 == 0)
                    epsilon *= 1.4142135;
                if (W[I, J] == 1e20)
                {
                    if (demgrid.Values[I, J] >= this.W[i, j] + epsilon)
                    {
                        this.W[I, J] = demgrid.Values[I, J];
                        Dry_upward_cell(I, J);
                    }
                }
            }
            depth--;

        }

        /// <summary>
        /// 经过本人改写后的Planchon算法
        /// </summary>
        /// <param name="DirectionMatrix"></param>
        /// <param name="dem"></param>
        public void CheckPit_Planchon_Modi(int[,] DirectionMatrix, GridLayer dem)
        {
            marks = new int[dem.rowCount, dem.colCount];
            marks2 = new int[demgrid.rowCount, demgrid.colCount];
            double[,] waterData=new double[dem.rowCount,dem.colCount];
            for (int row = 0; row < this.demgrid.rowCount; row++)
            {
                for (int col = 0; col < this.demgrid.colCount; col++)
                {
                    if (row == 0 || row == demgrid.rowCount - 1 || col == 0 || col == demgrid.colCount - 1)
                    {
                        waterData[row, col] = this.demgrid.Values[row, col];
                        marks[row, col] = 1;
                    }
                    else
                    {
                        waterData[row, col] = 1e20;
                        marks[row, col] = 0;
                    }
                }
            }
            for(int row=0;row<dem.rowCount;row++)
            {
                this.Process(waterData, row, 0);
                this.Process(waterData, row, dem.colCount - 1);
            }
            for(int col=0;col<dem.colCount;col++)
            {
                this.Process(waterData, 0, col);
                this.Process(waterData, dem.rowCount-1, col);
            }
            this.demgrid.Values = waterData;
        }
        private void Process(double[,] waterData, int row, int col)
        {
            Random random = new Random();
            ProcessCount++;
            Cell tCell = new Cell(row, col);
            Stack<Cell> stack = new Stack<Cell>();
            stack.Push(tCell);
            while (stack.Count > 0)
            {
                tCell = stack.Pop();
                int i = tCell.row;
                int j = tCell.col;
                //List<int> set = new List<int>();
                //for (int n = 0; n < 9; n++)
                //{
                //    set.Add(n);
                //}
                //for (int n = 0; n < 9; n++)
                //{
                //    int rnd=random.Next(set.Count);
                //    int index = set[rnd];
                //    set.RemoveAt(rnd);
                ////}
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    int I = i + index / 3 - 1;
                    int J = j + index % 3 - 1;
                    if (I < 0 || I >= this.demgrid.rowCount)
                        continue;
                    if (J < 0 || J >= this.demgrid.colCount)
                        continue;
                    double epsilon = 1e-10;
                    if (index % 2 == 0)
                        epsilon *= 1.4142135;

                    if (marks[I,J]==0 && marks2[I,J]!=ProcessCount)
                    {
                        stack.Push(new Cell(I, J));
                        marks2[I, J] = (int)ProcessCount;
                    }
                    if (this.demgrid.Values[i, j] >= waterData[I, J]+epsilon)
                    {
                        waterData[i, j] = this.demgrid.Values[i, j];
                        marks[i, j] = 1;
                        break;
                    }
                    if (waterData[i, j] > waterData[I, J]+epsilon)
                    {
                        waterData[i, j] = waterData[I, J]+epsilon;
                    }
                }
            }
            
        }
        private Cell NextCell(int R, int C, int i)
        {
            R += dR[i];
            C += dC[i];
            if (R < 0 || C < 0 || R >= this.demgrid.rowCount || C >= this.demgrid.colCount)
            {
                R += fR[i];
                C += fC[i];
            }
            if (R < 0 || C < 0 || R >= this.demgrid.rowCount || C >= this.demgrid.colCount)
            {
                return new Cell(-1,-1);
            }
            return new Cell(R, C);

        }
        private double depth = 0;
        private double[,] W;
        private bool something_done = false;
        
        int[] R0 = new int[8];
        int[] C0 = new int[8];
        int[] dR = new int[8];
        int[] dC = new int[8];
        int[] fR = new int[8];
        int[] fC = new int[8];
        public void Init_W()
        {            
            /*initial row*/
            int r= this.demgrid.rowCount - 2;
            int c = this.demgrid.colCount - 2;
            R0 = new int[] {1,r,1,r,1,r,1,r,1,r };           
            
            C0 = new int[] {1,c,c,1,c,1,1,c };            
            dR = new int[] {0,0,1,-1,0,0,1,-1 };
            dC = new int[] {1,-1,0,0,-1,1,0,0};

            fR = new int[] {1,-1,-r,r,1,-1,-r,r};
            
            fC = new int[] {-c,c,-1,1,c,-c,1,-1 };
           
            this.W=new double[demgrid.rowCount,demgrid.colCount];
            for (int i = 0; i < this.demgrid.rowCount; i++)
            {
                for (int j = 0; j < this.demgrid.colCount; j++)
                {
                    if (i == 0 || i == demgrid.rowCount - 1 || j == 0 || j == demgrid.colCount - 1)
                        this.W[i, j] = this.demgrid.Values[i, j];
                    else
                        this.W[i, j] = 1e20;
                }
            }
        }
           
        
        /// <summary>
        /// 计算水流累积矩阵，使用广度优先搜索
        /// </summary>
        /// <param name="MatrixDirection">已计算出的水流流向矩阵</param>
        /// <returns></returns>
        public int[,] FlowAccum_Old(int[,] DirectionLayer)
        {
            //ProgressBar bar = new ProgressBar();
            //bar.Show();
            //bar.Text = "正在计算累积流量矩阵";
            int[,] MatrixDirection = DirectionLayer;
            int[,] matrix = new int[demgrid.rowCount, demgrid.colCount];
            if (DirectionLayer == null)
            {
                MessageBox.Show("水流流向矩阵文件未读入，计算无法进行");
                return null;
            }
            
            for (int i = 0; i < demgrid.rowCount; i++)
            {
                for (int j = 0; j < demgrid.colCount; j++)
                {
                    //寻找当前像素(i,j)的所有流入像素

                    matrix[i,j]= this.AccumCountByPixel(MatrixDirection,i, j);    
                    //bar.progressBar1.Value=(int)(i*100.0/demgrid.rowCount);

                }
            }
            //bar.Close();
            return matrix;
        }
        public double[,] CalcuShade()
        {
            double[,] data = new double[this.demgrid.rowCount, demgrid.colCount];
            for (int i = 1; i < demgrid.rowCount-1; i++)
            {
                for (int j = 1; j < demgrid.colCount-1; j++)
                {
                    //计算向量
                    double x1 = 0;
                    double y1 = 1;
                    double z1 = demgrid.Values[i, j-1] - demgrid.Values[i - 1, j-1];

                    double x2 = 1;
                    double y2 = -1;
                    double z2 = demgrid.Values[i - 1, j] - demgrid.Values[i, j - 1];

                    double xx = y1 * z2 - z1 * y2;
                    double yy = z1 * x2 - x1 * z2;
                    double zz = x1 * y2 - y1 * x2;

                    double sum = Math.Sqrt(xx * xx + yy * yy + zz * zz);
                    if (sum < 0.000001)
                    {
                        xx = yy = zz = 1.0;
                    }
                    else
                    {
                        xx /= sum;
                        yy /= sum;
                        zz /= sum;
                    }
                    double[] normal = {xx,yy,zz };
                    double[] lightNormal = {1,2,1 };
                    data[i,j] = (int)(255 * (normal[0] * lightNormal[0] + normal[1] * lightNormal[1] + normal[2] * lightNormal[2]));
                }
            }
            return data;
        }
        private class XYComparer : IComparer<Cell>
        {
            public int Compare(Cell a,Cell b)
            {
                if (a.row > b.row)
                    return 1;
                if (a.row < b.row)
                    return -1;
                if (a.row == b.row)
                {
                    if (a.col > b.col)
                        return 1;
                    if (a.col < b.col)
                        return -1;
                }
                return 0; //a,b完全相等时为0
            }
            
        }
        
        /// <summary>
        /// 设定的方向编码方案如下：
        ///    i -------------
        ///    0 | 0 | 1 | 2 |
        ///      -------------
        ///    1 | 3 | 4 | 5 |
        ///      -------------
        ///    2 | 6 | 7 | 8 |
        ///      -------------
        ///        0   1   2 j
        /// 若当前像素(i,j)的流向编码为m，则此编码指向的相邻相素的行标和
        /// 列标分别为
        ///        I=i+m/3-1
        ///        J=j+m%3-1
        /// 当前像素某方向上的相邻像素的编码为a，则本像素位于相邻像素的方
        /// 向编码为8-a。其中编码4表示中心位置，在流向编码中应该避免使用
        /// 该码
        /// </summary>
        /// <param name="MatrixDirection"></param>
        /// <returns></returns>
        //public int[,] FlowAccumu_Liu()
        //{
        //    if (this.m_pDirection == null)
        //    {
        //        MessageBox.Show("未载入流向数据");
        //        return null;
        //    }

        //    ProgressBar bar = new ProgressBar();
        //    bar.Show();
        //    bar.Text = "正在生成流量累积矩阵";
        //    int row=demgrid.rowCount;
        //    int col=demgrid.colCount;
        //    int[,] DMatrix = Common.CopyMatrix(this.m_pDirection);   
            
        //    //存放累积值的矩阵
        //    int[,] MatrixAccumCounts = new int[demgrid.rowCount, demgrid.colCount];
        //    //把所有像素先放到一个set中，便于以后的快速搜索
        //    SortedDictionary<Cell, int> set = new SortedDictionary<Cell, int>();
        //    for (int i = 0; i < demgrid.rowCount; i++)
        //    {
        //        for (int j = 0; j < demgrid.colCount; j++)
        //        {
        //            set.Add(new Cell(i,j), -1);
        //        }
        //    }
            
        //    int count = 0;
        //    //int emptycount = 0;
        //    while (set.Count > 0)//当set不空时
        //    {                
        //        List<Cell> templist = new List<Cell>();
        //        //采用剥皮法统计累积流量
        //        foreach (Cell Key in set.Keys)
        //        {                   
        //            int i=Key.row;
        //            int j=Key.col;
        //            //如果该像元的记录无流入
        //            if (this.FindFlowInPixels(DMatrix, i, j).Count == 0)
        //            {
        //                //获取下游邻接像元，并对其累加流量
        //                Cell cell= FindFlowOutPixel(DMatrix, i, j);
        //                if(cell.row!=-1)
        //                    MatrixAccumCounts[cell.row, cell.col] += MatrixAccumCounts[i,j]+1;
        //                MatrixAccumCounts[i, j] += 1;
        //                DMatrix[i, j] = -1;//将其流向设为无
                        
        //                templist.Add(Key);
        //                count += 1;
        //                if (count % row == 0)
        //                {
        //                    bar.progressBar1.Value = (int)(count * 100.0 / row / col);
        //                    bar.Text="正在生成流量累积矩阵："+count.ToString()+" "+MatrixAccumCounts[i,j].ToString();
        //                }
        //            }
        //        }
        //        if (templist.Count == 0)
        //        {
        //            MessageBox.Show("可能因为数据的原因而分析失败！");
        //            break;
        //        }
        //        foreach (Cell mc in templist)
        //        {
        //            set.Remove(mc);
        //        }                
        //    }
        //    bar.Close();
        //    this.m_pAccumulation = MatrixAccumCounts;
        //    return MatrixAccumCounts;
        //}
        /// <summary>
        /// 递归函数，服务于 FlowAccumu_Liu5
        /// </summary>
        /// <returns></returns>
        private void FlowAccum(int row,int col)
        {
            if (this.m_pAccumulation[row, col] > 0) //如果面积为已知
            {
                return;
            }
            int[,] m = this.m_pDirection;//为了简化书写用m来代替MatrixDirection
            this.m_pAccumulation[row, col] = 1;
            //从八个方向搜索
            for (int index = 0; index < 9; index++)
            {
                if (index == 4)
                    continue;
                //获取index方向上的相邻像素的行标和列标
                int I = row + index / 3 - 1;
                int J = col + index % 3 - 1;

                if (I >= 0 && I < demgrid.rowCount && J >= 0 && J < demgrid.colCount && m[I, J] == 8 - index)  //如果相邻像素的位置指向本像素
                {
                    FlowAccum(I, J);
                    this.m_pAccumulation[row, col] += this.m_pAccumulation[I, J];
                }
            }
        }
        /// <summary>
        /// 计算累积流向，采用递归方法
        /// </summary>
        /// <returns></returns>
        public int[,] FlowAccumu_Liu5()
        {
            if (this.m_pDirection == null)
            {
                MessageBox.Show("未载入流向数据");
                return null;
            }

            //ProgressBar bar = new ProgressBar();
            //bar.Show();
            //bar.Text = "正在生成流量累积矩阵";
            int[,] DMatrix = MatrixFuncs.CopyMatrix(this.m_pDirection);

            //存放累积值的矩阵
            this.m_pAccumulation =new int[demgrid.rowCount, demgrid.colCount];
            for (int rowIndex = 0; rowIndex < demgrid.rowCount; rowIndex++)
            {
                this.FlowAccum(rowIndex, 0);//0列
                this.FlowAccum(rowIndex, demgrid.colCount - 1); //最后一列
            }
            for (int colIndex = 1; colIndex < demgrid.colCount - 1; colIndex++)
            {
                this.FlowAccum(0, colIndex);//0行
                this.FlowAccum(demgrid.rowCount - 1, colIndex);//最后一行
            }
            return this.m_pAccumulation;
        }
        public int[,] FlowAccumu_Liu3()
        {
            if (this.m_pDirection == null)
            {
                MessageBox.Show("未载入流向数据");
                return null;
            }

            //ProgressBar bar = new ProgressBar();
            //bar.Show();
            //bar.Text = "正在生成流量累积矩阵";
            int row = demgrid.rowCount;
            int col = demgrid.colCount;
            int[,] DMatrix = MatrixFuncs.CopyMatrix(this.m_pDirection);

            //存放累积值的矩阵
            int[,] MatrixAccumCounts = new int[demgrid.rowCount, demgrid.colCount];
            //把所有像素先放到一个set中，便于以后的快速搜索
            List<Cell> set = new List<Cell>();
            for (int i = 0; i < demgrid.rowCount; i++)
            {
                for (int j = 0; j < demgrid.colCount; j++)
                {
                    if(this.FindFlowInPixels(DMatrix,i,j).Count==0)
                        set.Add(new Cell(i, j));
                }
            }


            int count = 0;
            //int emptycount = 0;
            while (set.Count > 0)//当set不空时
            {
                List<Cell> NewSet = new List<Cell>();
                //List<Cell> templist = new List<Cell>();
                //采用剥皮法统计累积流量
                foreach (Cell Key in set)
                {
                    int i = Key.row;
                    int j = Key.col;
                    
                    //获取下游邻接像元，并对其累加流量
                    Cell cell = FindFlowOutPixel(DMatrix, i, j);
                    if (cell.row != -1)
                    {
                        MatrixAccumCounts[cell.row, cell.col] += MatrixAccumCounts[i, j] + 1;
                        
                        DMatrix[i, j] = -1;//将其流向设为无
                        if (this.FindFlowInPixels(DMatrix, cell.row, cell.col).Count == 0)
                            NewSet.Add(cell);                        
                    }
                    MatrixAccumCounts[i, j] += 1;
                    //count += 1;
                    //if (count % row == 0)
                    //{
                    //    bar.progressBar1.Value = (int)(count * 100.0 / row / col);
                    //    bar.Text = "正在生成流量累积矩阵：" + count.ToString() + " " + MatrixAccumCounts[i, j].ToString();
                    //}
                }
                set = NewSet;
            }
            //bar.Close();
            this.m_pAccumulation = MatrixAccumCounts;
            return MatrixAccumCounts;
        }
        public int[,] FlowAccumu_Liu2()
        {
            if (this.m_pDirection == null)
            {
                MessageBox.Show("未载入流向数据");
                return null;
            }

            //ProgressBar bar = new ProgressBar();
            //bar.Show();
            //bar.Text = "正在生成流量累积矩阵";
            int[,] DMatrix = MatrixFuncs.CopyMatrix(this.m_pDirection);

            //存放累积值的矩阵
            int[,] MatrixAccumCounts = new int[demgrid.rowCount, demgrid.colCount];

            for (int row = 0; row < demgrid.rowCount; row++)
            {
                for (int col = 0; col < demgrid.colCount; col++)
                {
                    int CurrentRow = row;
                    int CurrentCol = col;
                    while (true)
                    {
                        //获取下游像素的坐标
                        int Direction=DMatrix[CurrentRow,CurrentCol];
                        int NextRow = CurrentRow + Direction / 3 - 1;
                        int NextCol = CurrentCol + Direction % 3 - 1;
                        CurrentRow = NextRow;
                        CurrentCol = NextCol;
                        if (CurrentRow < 0 || CurrentRow > demgrid.rowCount - 1 || CurrentCol < 0 || CurrentCol > demgrid.colCount - 1)
                            break;
                        MatrixAccumCounts[CurrentRow, CurrentCol] += 1;
                    }
                    //bar.progressBar1.Value = (int)(row * 100.0 / demgrid.rowCount);
                }
            }
            this.m_pAccumulation = MatrixAccumCounts;
            //bar.Close();
            return MatrixAccumCounts;
        }

        public int[,] CalSubWaterShed(int[,] AccumMatrix)
        {
            
            if (this.m_pDirection == null)
            {
                MessageBox.Show("未载入流向数据");
                return null;
            }

            //ProgressBar bar = new ProgressBar();
            //bar.Show();
            //bar.Text = "正在生成流量累积矩阵";
            int row = demgrid.rowCount;
            int col = demgrid.colCount;
            int[,] DMatrix = MatrixFuncs.CopyMatrix(this.m_pDirection);

            //存放累积值的矩阵
            int[,] MatrixAccumCounts = new int[demgrid.rowCount, demgrid.colCount];
            //把所有像素先放到一个set中，便于以后的快速搜索
            List<Cell> set = new List<Cell>();
            for (int i = 0; i < demgrid.rowCount; i++)
            {
                for (int j = 0; j < demgrid.colCount; j++)
                {
                    if (this.FindFlowInPixels(DMatrix, i, j).Count == 0)
                        set.Add(new Cell(i, j));
                }
            }


            int count = 0;
            //int emptycount = 0;
            while (set.Count > 0)//当set不空时
            {
                List<Cell> NewSet = new List<Cell>();
                //List<Cell> templist = new List<Cell>();
                //采用剥皮法统计累积流量
                foreach (Cell Key in set)
                {
                    int i = Key.row;
                    int j = Key.col;

                    //获取下游邻接像元，并对其累加流量
                    Cell cell = FindFlowOutPixel(DMatrix, i, j);
                    if (cell.row != -1)
                    {
                        MatrixAccumCounts[cell.row, cell.col] += MatrixAccumCounts[i, j] + 1;

                        DMatrix[i, j] = -1;//将其流向设为无
                        if (this.FindFlowInPixels(DMatrix, cell.row, cell.col).Count == 0)
                            NewSet.Add(cell);
                    }
                    MatrixAccumCounts[i, j] += 1;
                    //count += 1;
                    //if (count % row == 0)
                    //{
                    //    bar.progressBar1.Value = (int)(count * 100.0 / row / col);
                    //    bar.Text = "正在生成流量累积矩阵：" + count.ToString() + " " + MatrixAccumCounts[i, j].ToString();
                    //}
                }
                set = NewSet;
            }
            //bar.Close();
            this.m_pAccumulation = MatrixAccumCounts;
            return MatrixAccumCounts;
        }
        /// <summary>
        /// 计算单个像素的累积流量
        /// </summary>
        /// <param name="MatrixDirection">流向矩阵</param>
        /// <param name="i">行标</param>
        /// <param name="j">列标</param>
        public int AccumCountByPixel(int[,] MatrixDirection, int i,int j)
        {
            int[,] m = MatrixDirection;//为了简化书写用m来代替MatrixDirection
            Queue<Cell> queue = new Queue<Cell>();//存放当前像素的所有上游流入像素
            this.PixelCount = 0;
            queue.Enqueue(new Cell(i, j));
            
            while (queue.Count > 0)
            {
                Cell cell = queue.Dequeue();
                i = cell.row;
                j = cell.col;
                //从八个方向搜索
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    //获取index方向上的相邻像素的行标和列标
                    int I = i + index / 3 - 1;
                    int J = j + index % 3 - 1;

                    if (I >= 0 && I < demgrid.rowCount && J >= 0 && J < demgrid.colCount && m[I, J] == 8 - index)  //如果相邻像素的位置指向本像素
                    {
                        queue.Enqueue(new Cell(I,J));
                        this.PixelCount += 1;
                    }
                }
               
            }
            return this.PixelCount;
        }
        /// <summary>
        /// 从整数矩阵中查找最大的元素值
        /// </summary>
        /// <param name="MatrixAccum"></param>
        /// <returns></returns>
        public int FindMaxValue(int[,] Matrix)
        {
            int max = -99999;
            for (int i = 0; i < this.demgrid.rowCount; i++)
            {
                for (int j = 0; j < this.demgrid.colCount; j++)
                {
                    if (Matrix[i, j] > max)
                        max = Matrix[i, j];
                }
            }
            return max;
        }
        /// <summary>
        /// 生成流向位图
        /// </summary>
        /// <param name="matrix"></param>
        /// <returns></returns>
        public Bitmap CreateFlowDirectImage(int[,] DirectionLayer)
        {
            ProgressBar bar = new ProgressBar();
            bar.Show();
            bar.Text = "正在生成位图";
            int[,] matrix = DirectionLayer;
            Bitmap bitmap = new Bitmap(demgrid.colCount, demgrid.rowCount);
            for (int i = 0; i < demgrid.rowCount; i++)
            {
                for (int j = 0; j < demgrid.colCount; j++)
                {
                    if (matrix[i, j] > 0)
                    {
                        int v = matrix[i, j];
                        if (v < 7)
                            v += 8;
                        v = (v - 7) * 32;
                        bitmap.SetPixel(j, i, Color.FromArgb(v, v, v));
                        bar.progressBar1.Value = (int)(i * 100.0 / demgrid.rowCount);
                    }
                    else
                        bitmap.SetPixel(j, i, Color.Black);
                }
            }
            bar.Close();
            return bitmap;
        }
        /// <summary>
        /// 通过累积流量矩阵生成一个水系位图
        /// </summary>
        /// <param name="MatrixAccum">累积流量矩阵</param>
        /// <param name="threshold">生成水系时用的累积流量门限值</param>
        /// <returns></returns>
        public Bitmap CeateRiverSysImage(int[,] MatrixAccum, int threshold)
        {
            int max=this.FindMaxValue(MatrixAccum);
            Bitmap bitmap = new Bitmap(demgrid.colCount, demgrid.rowCount);
            for (int i = 0; i < demgrid.rowCount; i++)
            {
                for (int j = 0; j < demgrid.colCount; j++)
                {
                    if (MatrixAccum[i, j] > threshold)
                        bitmap.SetPixel(j, i, Color.White);
                    else
                        bitmap.SetPixel(j, i, Color.Black);
                }
            }
            return bitmap;
        }
        /// <summary>
        /// 寻找当前像素的下游流入像素
        /// </summary>
        /// <param name="MatrixDirection">流向矩阵</param>
        /// <param name="i">当前像素的行号</param>
        /// <param name="j">当前像素的列号</param>
        /// <returns></returns>
        private Cell FindFlowOutPixel(int[,] MatrixDirection, int i,int j)
        {
            int I = i + MatrixDirection[i, j] / 3 - 1;
            if (I >= demgrid.rowCount || I < 0)  //下标越界
                return new Cell(-1,-1);
            int J = j + MatrixDirection[i, j] % 3 - 1;
            if (J >= demgrid.colCount || J < 0)//下标越界
                return new Cell(-1,-1);
            return new Cell(I, J);
        }
        /// <summary>
        /// 计算流入当前像素(x,y)的相邻像素数
        /// </summary>
        /// <param name="MatrixDirection">流向矩阵</param>
        /// <param name="i">行号，注意它只对应图像中的纵坐标</param>
        /// <param name="j">列号，注意它只对应图像中的横坐标</param>
        /// <returns></returns>
        private List<Cell> FindFlowInPixels(int[,] MatrixDirection, int i,int j)
        {
            int[,] m=MatrixDirection;
            List<Cell> cells = new List<Cell>();
            //遍历8个方向
            for (int index = 0; index < 9; index++)
            {
                if (index == 4)
                    continue;
                //获取index方向上的相邻像素的行标和列标
                int I = i + index / 3 - 1;
                int J = j + index % 3 - 1;
                
                if (I>=0 && I<demgrid.rowCount && J>=0 && J<demgrid.colCount)
                {
                    int direct = m[I, J];
                    if(direct == 8 - index)  //如果相邻像素的位置指向本像素
                        cells.Add(new Cell(I, J));
                }
            }
            return cells;
        }
        /// <summary>
        /// 获取当前像素上游那些超过域值threshold的所有流入像素
        /// </summary>
        /// <param name="MatrixDirection">流向值网格</param>
        /// <param name="AccumLayer">累积流量网格</param>
        /// <param name="i">当前像素行号</param>
        /// <param name="j">当前像素列号</param>
        /// <param name="threshold">累积流量门限值</param>
        /// <returns></returns>
        private List<Cell> FindFlowInPixels(int[,] MatrixDirection, int[,] AccumLayer, int i, int j,int threshold)
        {
            int[,] m = MatrixDirection;
            int[,] am=AccumLayer;
            List<Cell> cells = new List<Cell>();
            //遍历8个方向
            for (int index = 0; index < 9; index++)
            {
                if (index == 4)
                    continue;
                //获取index方向上的相邻像素的行标和列标
                int I = i + index / 3 - 1;
                int J = j + index % 3 - 1;

                if (I >= 0 && I < demgrid.rowCount && J >= 0 && J < demgrid.colCount && m[I, J] == 8 - index)  //如果相邻像素的位置指向本像素
                {
                    if(am[I,J]>threshold)
                        cells.Add(new Cell(I, J));
                }
            }
            return cells;
        }
        /// <summary>
        /// Strahler河流分级
        /// </summary>
        /// <param name="cellX"></param>
        /// <param name="cellY"></param>
        /// <param name="NodeArray"></param>
        /// <returns></returns>
        public int[,] FlowStrahlerOrder(int threshold)
        {
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在生成河网";
            bar.Show();
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;

            int[,] MMM = new int[rows, cols];

            //用于存放属于河流的像素，及像素上的级别
            RiverNet = new SortedDictionary<Cell, int>();
            int[,] matrix =this.m_pAccumulation;

            //获取河流像素
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    if (matrix[i, j] > threshold)
                        RiverNet.Add(new Cell(i, j), matrix[i, j]);
                }
            }

            
            //从第0级开始搜索
            //int index = 0;
            while (true)
            {
                List<Cell> cells = new List<Cell>();
                foreach (KeyValuePair<Cell, int> pair in RiverNet)
                {
                    Cell cell = pair.Key;
                    if (MMM[cell.row, cell.col] != 0)//如果该像素已有级别，则跳过不管
                        continue;

                    List<Cell> PreCells = this.FindFlowInPixels(this.m_pDirection,this.m_pAccumulation, pair.Key.row, pair.Key.col, threshold);
                    //求出流入像素中还未有级别的像素个数

                    if (PreCells.Count == 0) //如果没有流入像素，则更新级别值为1
                    {
                        MMM[cell.row, cell.col] = 1;
                        cells.Add(cell);
                        continue;
                    }
                    //如果有流入像素，先求流入像素中未更新级别值的个数                    
                    int mcount = 0;  //流入像素中级别值未更新的像素个数
                    int maxvalue = -1; //流入像素中最大级别值
                    int maxvalueCount = 0; //具有最大级别值的像素的个数
                    foreach (Cell mc in PreCells)
                    {
                        int temp = MMM[mc.row, mc.col];
                        if (temp == 0)
                            mcount += 1;
                        if (maxvalue < temp)
                        {
                            maxvalue = temp;
                            maxvalueCount = 1;
                        }
                        else if (maxvalue == temp)
                            maxvalueCount += 1;
                    }

                    if (mcount == 0) //如果全部更新完毕，就更新当前像素值
                    {
                        if (PreCells.Count == 1)  //如果只有一个流入像素，则级别与流入相同
                        {
                            MMM[cell.row, cell.col] = MMM[PreCells[0].row, PreCells[0].col];
                            cells.Add(cell);
                        }
                        else  //级别值有多个时
                        {
                            if (maxvalueCount > 1)//有至少两个相同最高级别的流入像素
                            {
                                MMM[cell.row, cell.col] = maxvalue + 1;
                                cells.Add(cell);
                            }
                            else
                            {
                                MMM[cell.row, cell.col] = maxvalue;
                                cells.Add(cell);
                            }
                        }
                    }
                }
                if (cells.Count == 0)
                    break;
                foreach (Cell cc in cells)  //删去
                {
                    RiverNet.Remove(cc);
                }

            }

            bar.Close();
            this.m_pStrahlerOrd = MMM;
            return MMM;
        }
        /// <summary>
        /// 提取河道
        /// </summary>
        /// <param name="AccumuMatrix">流量累积矩阵</param>
        /// <param name="FlowDirectionLayer">流向矩阵</param>
        /// <param name="threshold"></param>
        /// <returns></returns>
        public int[,] FlowStrahlerOrder1(int[,] AccumuMatrix,int threshold)
        {
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;

            int[,] MMM = new int[rows, cols];

            //用于存放属于河流的像素，及像素上的级别
            RiverNet = new SortedDictionary<Cell, int>();
            int[,] matrix = AccumuMatrix;

            //获取河流像素
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    if (matrix[i, j] > threshold)
                        RiverNet.Add(new Cell(i, j), matrix[i, j]);
                }
            }
                       
            return MMM;

        }

        public bool FlowMaximumLength(int cellX, int cellY)
        {
            return true;
        }
        /*======================================================================+
    +																		+
    +	函数名：	FlowTotalLength											+
    +																		+
    +	输入：		cellX - 河道出口行号；cellY - 河道出口列号；			+
    +																		+
    +	输出：		true - 分析完成；false - 分析失败；						+
    +				河道总长度矩阵 - m_pTotalFlowLength；					+
    +																		+
    +	所需矩阵：	河网水系矩阵 - m_pStreamNet；							+
    +				水流流向矩阵 - m_pDirection;							+
    +																		+
    +	功能：		提取河道总长度(从河道末端栅格计算到流域出口栅格)		+
    +																		+
    +	原理：		对已经按阈值提取的河道计算河道长度划分,所以河道			+
    +				依赖于最小集水面积阈值的设置。							+
    +																		+
    +======================================================================*/
        bool FlowTotalLength(int cellX, int cellY)
        {
            //string title1, title2;
            //if (m_pDirection == null)
            //{
            //    //title1="水流流向数据文件*_fld";
            //    //title2="读入水流流向数据文件";
            //    //if(!DirectRead(title1,title2,"asc",m_pDirection))	
            //    //{
            //    //    ErrorReport(1);
            //    //    return false;
            //    //}
            //}
            //if (m_pStreamNet == null)
            //{
            //    //title1="河网水系矩阵文件*_snk";
            //    //title2="读入河网水系矩阵文件";
            //    //if(!DirectRead(title1,title2,"asc",&m_pStreamNet))	
            //    //{
            //    //    ErrorReport(6);
            //    //    return false;
            //    //}
            //}



            //int PixelsCount = m_row * m_col;
            //int[] ii = new int[PixelsCount];
            //int[] jj = new int[PixelsCount];
            //for (int i = 0; i < PixelsCount; i++)
            //{
            //    ii[i] = 0;
            //    jj[i] = 0;
            //}

            //ii[0] = cellX;
            //jj[0] = cellY;

            //int inum = 0;
            //int totin = 1;
            //while (inum < totin)
            //{
            //    //totin += PixelFlowIn(ii[inum], jj[inum], ii, jj, totin, m_pTotalFlowLength);
            //    inum++;
            //}

            //int[]i8 = new int[8];
            //int[]j8 = new int[8];
            //int iInGrid = 0;
            //int pos;

            //for (int k = inum - 1; k >= 0; k--)
            //{
            //    if (m_pStreamNet[ii[k]][jj[k]] != 1.0)
            //        continue;
            //    pos = 0;
            //    iInGrid = PixelFlowIn(ii[k], jj[k], i8, j8, pos, m_pStreamNet, 1.0);
            //    if (iInGrid == 0)
            //    {
            //        m_pTotalFlowLength[ii[k]][jj[k]] = 0.0;
            //    }
            //    else if (iInGrid == 1)
            //    {
            //        m_pTotalFlowLength[ii[k]][jj[k]] = m_pTotalFlowLength[i8[0]][j8[0]]
            //            + GridFlowLength((int)m_pDirection[i8[0]][j8[0]], demgrid.resolusion);
            //    }
            //    else if (iInGrid > 1)
            //    {
            //        for (int j = 0; j < iInGrid; j++)
            //        {
            //            m_pTotalFlowLength[ii[k]][jj[k]] = m_pTotalFlowLength[ii[k]][jj[k]]
            //                + m_pTotalFlowLength[i8[j]][j8[j]] +
            //                GridFlowLength((int)m_pDirection[i8[j]][j8[j]], demgrid.resolusion);
            //        }
            //    }
            //}

            //if (ii) { delete[] ii; ii = null; }
            //if (jj) { delete[] jj; jj = null; }
            //if (i8) { delete[] i8; i8 = null; }
            //if (j8) { delete[] j8; j8 = null; }
            return true;
        }


        

        /*======================================================================+
        +																		+
        +						 子流域边界提取计算模块							+
        +																		+
        +				      （已经通过测试，2005-03-30）						+
        +				   （计算原理见 模型系统设计书 第二章 2.4.6）			+
        +																		+
        +======================================================================*/
    //    public bool WaterShed(int cellX, int cellY)
    //    {
    //        if(m_pDirection==null)		
    //{
    //    string title1,title2;
    //    title1="水流流向数据文件*_fld";
    //    title2="读入水流流向数据文件";
    //    if(!DirectRead(title1,title2,"asc",m_pDirection))	
    //    {
    //        ErrorReport(1);
    //        return false;
    //    }
    //}
    //if(m_pWatershed==null)
    //    m_pWatershed=new double[this.demgrid.row];
 	
	
    //int totnum = m_row * m_col;		
    //int *ii = new int[totnum];
    //int *jj = new int[totnum];
    //for(int i=0;i<totnum;i++)
    //{
    //    ii[i]=0;
    //    jj[i]=0;
    //}

    //ii[0] = cellX;
    //jj[0] = cellY;
	
    //int inum=0;
    //int totin = 1;
    //while (inum<totin) 
    //{
    //    totin += PixelFlowIn(ii[inum],jj[inum],ii,jj,totin,m_pWatershed);
    //    inum++;		
    //}
	
    //for(int k = inum - 1; k >= 0 ; k--)
    //{
    //    m_pWatershed[ii[k]][jj[k]] = 1;				
    //}	

    //if(ii){delete[]ii;ii=null;}
    //if(jj){delete[]jj;jj=null;}
    //return true;
    //    }

        
        
        /// <summary>
        /// 计算水流流向矩阵(及其反向矩阵)
        /// </summary>
        /// <param name="bStat"></param>
        /// <returns></returns>
        public int[,] GridGradientDirection()
        {
            if(this.demgrid==null || this.demgrid.Values==null)
            {
                MessageBox.Show("未载入DEM数据");
                return null;
            }
            
            ProgressBar bar = new ProgressBar();
            bar.Show();
            bar.Text = "正在计算流向矩阵";
            int[,] DirectiontMatrix = new int[demgrid.rowCount, demgrid.colCount];
            for (int i = 1; i < demgrid.rowCount-1; i++)
            {
                for (int j = 1; j <demgrid.colCount-1; j++)
                {
                    if (demgrid.Values[i, j] < 0)
                    {
                        DirectiontMatrix[i, j] = -2;
                    }
                    //得到像元i,j的流向值
                    else
                    {
                        int value = PixelFlowOut(this.demgrid.Values, i, j);
                        DirectiontMatrix[i, j] = value;
                    }
                }
                bar.progressBar1.Value = (int)(i * 100.0 / demgrid.rowCount);
            }
            for (int j = 0; j < demgrid.colCount; j++)//顶行及底行像元
            {
                DirectiontMatrix[0, j] = 1;
                DirectiontMatrix[demgrid.rowCount - 1, j] = 7;
            }


            for (int i = 0; i < demgrid.rowCount; i++)//左边及右边像元
            {
                DirectiontMatrix[i, 0] = 3;
                DirectiontMatrix[i, demgrid.colCount - 1] = 5;
            }
                           
            
            bar.Close();
            this.m_pDirection = DirectiontMatrix;
            return DirectiontMatrix;
        }

        /// <summary>
        /// 当流向较多时需要处理为1个流向
        /// </summary>
        /// <param name="i8">数组，存储八个方向的流向标志</param>
        /// <returns>返回最后决定的方向</returns>
        private int GridMultiDirection(List<int> MaxIndexValues)
        {
            if (MaxIndexValues.Count == 1)
                throw new Exception("");
            if (MaxIndexValues.Count == 2)
            {
                Random random = new Random(123456);
                int index = (int)(random.NextDouble() * 2);
                return MaxIndexValues[index];
            }
            if (MaxIndexValues.Contains(0) && MaxIndexValues.Contains(1) && MaxIndexValues.Contains(2))
                return 1;
            if (MaxIndexValues.Contains(1) && MaxIndexValues.Contains(2) && MaxIndexValues.Contains(5))
                return 2;
            if (MaxIndexValues.Contains(0) && MaxIndexValues.Contains(3) && MaxIndexValues.Contains(6))
                return 3;
            if (MaxIndexValues.Contains(2) && MaxIndexValues.Contains(5) && MaxIndexValues.Contains(8))
                return 5;
            if (MaxIndexValues.Contains(3) && MaxIndexValues.Contains(6) && MaxIndexValues.Contains(7))
                return 6;
            if (MaxIndexValues.Contains(6) && MaxIndexValues.Contains(7) && MaxIndexValues.Contains(8))
                return 7;
            if (MaxIndexValues.Contains(7) && MaxIndexValues.Contains(8) && MaxIndexValues.Contains(5))
                return 8;
            return -1;
        }

        
        /// <summary>
        /// 获取DEM中坐标为cellX和cellY的象元上的流向
        /// </summary>
        /// <param name="v">高程矩阵</param>
        /// <param name="cellX">行标</param>
        /// <param name="cellY">列标</param>
        /// <returns>返回流向编码</returns>    
        public int PixelFlowOut(double[,] v, int cellX, int cellY)
        {            
            double CurrentValue = v[cellX, cellY];
            if (CurrentValue<1e-6)
                return -2;
            //下面寻找当前象素周围最大的坡度值，或者称梯度值
            
            int MaxValueIndex = -1;
            double dMax = 0;
            for (int index = 0; index < 9; index++)
            {
                if (index == 4)
                    continue;
                int I = cellX + index / 3 - 1;
                int J = cellY + index % 3 - 1;
                if (I < 0 || I >= demgrid.rowCount)
                    continue;
                if (J < 0 || J >= demgrid.colCount)
                    continue;

                double grad = 0;
                if (index % 2 == 0) //偶数为对角线方向编码
                    grad = (CurrentValue - v[I, J]) / 1.41421;
                else
                    grad = CurrentValue - v[I, J];
                if (grad > dMax)
                {
                    dMax = grad;
                    MaxValueIndex=index;
                }
                if(marks2!=null && grad==0 && marks2[I,J]<marks2[cellX,cellY])
                {
                    MaxValueIndex=index;
                }
            }
            return MaxValueIndex;
        }
       

        int PixelRouteOrder(int cellX, int cellY, int ii, int jj, int pos, double pValue)
        {
            return 0;
        }

        /*======================================================+
        +						其它函数定义					+
        +======================================================*/

        double GetGridValue(double pValue, int x, int y)
        {
            return 0;
        }
        void GetMinMaxValue(double pValue, double dMax, double dMin)
        {
        }
        

        double GridFlowLength(int nDirect, double resolusion)
        {
            double dret = 0.0;
            //if (nDirect == d8[0] || nDirect == d8[2] || nDirect == d8[4] || nDirect == d8[6])
            //    dret = resolusion;
            //if (nDirect == d8[1] || nDirect == d8[3] || nDirect == d8[5] || nDirect == d8[7])
            //    dret = resolusion *Math.Sqrt(2);
            return dret;
        }


        /// <summary>
        /// 计算流程长度矩阵(用于滞时计算)，同时生成汇流次序数据表。(汇流次序等级矩阵被流程长度矩阵代替)
        /// </summary>
        /// <param name="MatrixDirection"></param>
        /// <param name="AccumuMatrix"></param>
        /// <returns>流程长度矩阵</returns>
        public double[,] FlowCalcuOrder()
        {
            int[,] AccumuMatrix = this.m_pAccumulation;
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流次序等级矩阵";
            bar.Show();
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;
            int MaxRow = -1;
            int MaxCol = -1;
            int MaxValue = -1;
            //从边界上查找累积流量最大的像素，它就是河流的出口
            for (int row = 0; row < rows; row++)
            {
                if (AccumuMatrix[row, 0] > MaxValue)
                {
                    MaxValue = AccumuMatrix[row, 0];
                    MaxRow = row;
                    MaxCol = 0;
                }
                if (AccumuMatrix[row, cols - 1] > MaxValue)
                {
                    MaxValue = AccumuMatrix[row, cols - 1];
                    MaxRow = row;
                    MaxCol = cols - 1;
                }
            }

            for (int col = 0; col < cols; col++)
            {
                if (AccumuMatrix[0, col] > MaxValue)
                {
                    MaxValue = AccumuMatrix[0, col];
                    MaxRow = 0;
                    MaxCol = col;
                }
                if (AccumuMatrix[rows - 1, col] > MaxValue)
                {
                    MaxValue = AccumuMatrix[rows - 1, col];
                    MaxRow = rows - 1;
                    MaxCol = col;
                }
            }
            //记录高程最低的像素
            FindedCell = new Cell(MaxRow, MaxCol);
            //下面从流域中高程最低像素开始计算汇流次序
            double[,] result = new double[rows, cols];

            Queue<Cell> queue = new Queue<Cell>();
            queue.Enqueue(new Cell(MaxRow, MaxCol));
            this.CalcuOrderList = new List<Cell>();
            this.CalcuOrderList.Add(new Cell(MaxRow, MaxCol));
            int CountCells = rows * cols;
            int Count = 0;
            while (queue.Count > 0)
            {
                Cell CurCell = queue.Dequeue();
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    int R = CurCell.row + index / 3 - 1;
                    int C = CurCell.col + index % 3 - 1;
                    if (R < 0 || R >= this.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= this.demgrid.colCount)
                        continue;
                    if (this.m_pDirection[R, C] == 8 - index)
                    {
                        if (index % 2 == 0) //斜对角方向
                            result[R, C] = result[CurCell.row, CurCell.col] + 1.41421356;
                        else //正方向
                            result[R, C] = result[CurCell.row, CurCell.col] + 1;
                        queue.Enqueue(new Cell(R, C));
                        CalcuOrderList.Add(new Cell(R, C));
                        Count++;
                        if (Count % 100 == 0)
                        {
                            bar.progressBar1.Value = (int)(100.0 * Count / CountCells);
                        }
                    }
                }
            }
            bar.Close();
            CalcuOrderList.Reverse();
            this.m_pGridRouteLength = result;
            return result;

        }
        /// <summary>
        /// 指定河流出口断面后计算流程长度矩阵(用于滞时计算)，同时生成汇流次序数据表。
        /// </summary>
        /// <param name="row">出口断面所在栅格行号</param>
        /// <param name="col">出口断面所在栅格列号</param>
        /// <returns>流程长度矩阵(路径长度单位是栅格个数)</returns>
        public double[,] FlowCalcuOrder2(int row,int col)
        {
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流次序等级矩阵";
            bar.Show();
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;
            
            //记录高程最低的像素
            FindedCell = new Cell(row,col);
            //下面从流域中高程最低像素开始计算汇流次序
            double[,] RouteLength = new double[rows, cols];
            this.m_pRouteOrd = new double[rows, cols];
            Queue<Cell> queue = new Queue<Cell>();
            queue.Enqueue(FindedCell);
            RouteLength[FindedCell.row,FindedCell.col] = 1;
            this.CalcuOrderList = new List<Cell>();
            
            
            this.CalcuOrderList.Add(FindedCell);
            int CountCells = rows * cols;
            int Count = 0;
            while (queue.Count > 0)
            {
                Cell CurCell = queue.Dequeue();
                
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    int R = CurCell.row + index / 3 - 1;
                    int C = CurCell.col + index % 3 - 1;
                    if (R < 0 || R >= this.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= this.demgrid.colCount)
                        continue;
                    if (this.m_pDirection[R, C] == 8 - index)
                    {
                        if (index % 2 == 0) //斜对角方向
                            RouteLength[R, C] = RouteLength[CurCell.row, CurCell.col] + 1.41421356;
                        else //正方向
                            RouteLength[R, C] = RouteLength[CurCell.row, CurCell.col] + 1;
                        queue.Enqueue(new Cell(R, C));
                        CalcuOrderList.Add(new Cell(R, C));
                        this.m_pRouteOrd[R, C] = CalcuOrderList.Count;
                        Count++;
                        if (Count % 100 == 0)
                        {
                            bar.progressBar1.Value = (int)(100.0 * Count / CountCells);
                        }
                    }
                }
            }
            bar.Close();            
            CalcuOrderList.Reverse();
            //this.SaveCalcuOrder(HydroSimulate.path + "\\SysPara\\CalcuOrder.txt");
            //Common.SaveMatrix(RouteLength, "C:\\routelength.txt");
            this.m_pGridRouteLength = RouteLength;
            return RouteLength;

        }
        
        /// <summary>
        /// 仅计算汇流次序表(2016.12.4)
        /// </summary>
        /// <param name="row">出口断面所在栅格行号</param>
        /// <param name="col">出口断面所在栅格列号</param>
        /// <returns>流程长度矩阵(路径长度单位是栅格个数)</returns>
        public List<double[]> FlowCalcuOrder3(int row,int col)
        {
            
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;
            
            //记录高程最低的像素
            FindedCell = new Cell(row,col);
            //下面从流域中高程最低像素开始计算汇流次序
           
            
            Queue<Cell> queue = new Queue<Cell>();
            queue.Enqueue(FindedCell);
            List<double[]> Orders=new List<double[]>();
            Orders.Add(new double[]{row,col,-1,-1,0});
            double dist=1;
            
            while (queue.Count > 0)
            {
                Cell CurCell = queue.Dequeue();
                
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    int R = CurCell.row + index / 3 - 1;
                    int C = CurCell.col + index % 3 - 1;
                    if (R < 0 || R >= this.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= this.demgrid.colCount)
                        continue;
                    if (this.m_pDirection[R, C] == 8 - index)
                    {
                        if (index % 2 == 0) //斜对角方向
                            dist= 1.41421356;
                        else //正方向
                            dist= 1;
                        queue.Enqueue(new Cell(R, C));
                        Orders.Add(new double[]{R,C,CurCell.row,CurCell.col,dist});                        
                    }
                }
            }
           
            Orders.Reverse();
            StreamWriter sw=new StreamWriter("C:/FlowOrder.txt");
            for(int i=0;i<Orders.Count;i++)
            {
            	double[] values=Orders[i];
            	for(int j=0;j<values.Length;j++)
            	{
            		sw.Write(values[j]+"\t");
            	}
            	sw.WriteLine();
            }
            sw.Close();
            //this.SaveCalcuOrder(HydroSimulate.path + "\\SysPara\\CalcuOrder.txt");
            //Common.SaveMatrix(RouteLength, "C:\\routelength.txt
            return Orders;

        }
        /// <summary>
        /// 划分子流域
        /// </summary>
        /// <param name="row">出口断面所在栅格行号</param>
        /// <param name="col">出口断面所在栅格列号</param>
        /// <returns>流程长度矩阵(路径长度单位是栅格个数)</returns>
        public int[,] SubWaterShed(int row, int col,double threshold)
        {
            if (this.m_pAccumulation == null)
                throw new Exception("还没有生成累积流量矩阵");
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流次序等级矩阵";
            bar.Show();
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;
            int CurrentID = 0;
            //记录高程最低的像素
            FindedCell = new Cell(row, col);
            //下面从流域中高程最低像素开始计算汇流次序
            int[,] SubWaterSheds = new int[rows, cols];
            //this.m_pRouteOrd = new double[rows, cols];
            Queue<Cell> queue = new Queue<Cell>();
            queue.Enqueue(FindedCell);
            SubWaterSheds[FindedCell.row, FindedCell.col] = CurrentID;
            
            //int CountCells = rows * cols;
            int Count = 0;
            int SubCount = 1;
            while (queue.Count > 0)
            {
                Cell CurCell = queue.Dequeue();
                int flowInCount=this.FindFlowInPixels(this.m_pDirection, CurCell.row, CurCell.col).Count;
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;
                    int R = CurCell.row + index / 3 - 1;
                    int C = CurCell.col + index % 3 - 1;
                    if (R < 0 || R >= this.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= this.demgrid.colCount)
                        continue;
                    if (this.m_pDirection[R, C] == 8 - index)
                    {                        
                        queue.Enqueue(new Cell(R, C));
                        //如果(R,C)的累积流量大于阈值，且与Curcell累积流量的差大于等于阈值，则(R,C)是一个新子流域,
                        //否则为老子流域
                        double thr0=m_pAccumulation[R,C];
                        double thr1=m_pAccumulation[CurCell.row,CurCell.col];
                        if (thr0 >= threshold && thr1-thr0>=threshold)
                        {
                            SubWaterSheds[R, C] = SubWaterSheds[CurCell.row, CurCell.col] + 1;
                            SubCount += 1;
                        }
                        else
                        {
                            SubWaterSheds[R, C] = SubWaterSheds[CurCell.row, CurCell.col];
                        }
                        //Count++;
                        //if (Count % 100 == 0)
                        //{
                        //    bar.progressBar1.Value = (int)(100.0 * Count / CountCells);
                        //}
                    }
                }
            }
            bar.Close();
            this.m_pSubWatershed = SubWaterSheds;
            return SubWaterSheds;

        }
        public void SaveCalcuOrder(string FileName)
        {
            if (this.CalcuOrderList == null)
            {
                MessageBox.Show("你还没有生成汇流次序表");
                return;
            }
            StreamWriter sw = new StreamWriter(FileName);
            foreach (Cell cell in this.CalcuOrderList)
            {
                sw.Write(cell.row + "\t");
                sw.WriteLine(cell.col + "\t");
            }
            sw.Close();
        }
        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +															+
        +					栅格汇流时间计算						+
        +															+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +					  L          L							+
        +				tc = --- = ------------						+
        +					  V		dKv*Sqrt(S)						+
        +															+
        +				L   -- 栅格的水流路径长度；					+
        +				S   -- 栅格的平均坡度；						+
        +				dKV -- 速度常数，包含糙率、水力半径等因素	+
        +						对水流的影响(外部输入，需率定)		+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
        /// <summary>
        /// 计算栅格汇流时间
        /// </summary>
        /// <param name="SurfKv">地表径流速度常数</param>
        /// <param name="LatKv">壤中流速度常数</param>
        /// <param name="BaseKv">基流速度常数</param>
        /// <param name="RiverKv">河道径流速度常数</param>
        public void CalcRouteTime(double SurfKv,double LatKv,double BaseKv,double RiverKv)
        {
            
            int rows = demgrid.rowCount;
            int cols = demgrid.colCount;
            this.m_pLatQRouteTime = new double[rows, cols];
            this.m_pSurfQRouteTime = new double[rows, cols];
            this.m_pBaseQRouteTime = new double[rows, cols];
            double LowestValue = this.demgrid.Values[this.FindedCell.row, this.FindedCell.col];
            double[,] DEM=this.demgrid.Values;
            for (int row = 0; row < rows; row++)
            {
                for (int col = 0; col < cols; col++)
                {
                    //计算当前像素与最低像素的高程差
                    double heightDiff = DEM[row, col] - LowestValue;
                    //平均坡度
                    double length=this.m_pGridRouteLength[row,col]*this.demgrid.resolution;
                    if (length == 0)
                        continue;
                    if (heightDiff == 0)//指流域出口处
                    {
                        this.m_pSurfQRouteTime[row, col] = 0;
                        this.m_pLatQRouteTime[row, col] = 0;
                        this.m_pBaseQRouteTime[row, col] = 0;
                    }
                    else
                    {
                        double avgSlope = heightDiff / length;
                        double SqrtAvgSlope = Math.Sqrt(avgSlope);

                        //如果坡度太小，则等于不能汇流，用-1来表示
                        if (SqrtAvgSlope < 1e-20)
                        {
                            this.m_pSurfQRouteTime[row, col] = -1;
                            this.m_pLatQRouteTime[row, col] = -1;
                            this.m_pBaseQRouteTime[row, col] = -1;
                        }
                        else
                        {
                            
                            this.m_pSurfQRouteTime[row, col] = length / SqrtAvgSlope / SurfKv/3;
                            this.m_pLatQRouteTime[row, col] = length / SqrtAvgSlope / LatKv/3;
                            this.m_pBaseQRouteTime[row, col] = length / SqrtAvgSlope / BaseKv/3;
                        }
                    }
                }
            }
            MatrixFuncs.SaveMatrix(this.m_pSurfQRouteTime, "C://surftime.txt");
            MatrixFuncs.SaveMatrix(this.m_pLatQRouteTime, "C://lattime.txt");
            MatrixFuncs.SaveMatrix(this.m_pBaseQRouteTime, "C://basetime.txt");
            
        }
        
        public GridLayer demgrid;
        /// <summary>
        /// 汇流计算次序列表
        /// </summary>
        public List<Cell> CalcuOrderList;
        int m_col, m_row;
        double m_nodata;

        //DEM河道参数(自动提取)
        /// <summary>
        /// 坡度矩阵
        /// </summary>
      public  double[,] m_pSlope;				//坡度
        /// <summary>
        /// 坡向距阵
        /// </summary>
       public double[,] m_pAspect;			//坡向
        /// <summary>
        /// 洼地矩阵
        /// </summary>
       public double[,] m_pDepression;		//洼地
        /// <summary>
        /// 水流流向矩阵
        /// </summary>
        public int[,] m_pDirection;			//水流流向矩阵
        /// <summary>
        /// 水流累积量矩阵
        /// </summary>
       public int[,] m_pAccumulation;		//水流累积量矩阵
        /// <summary>
        /// 河网水系矩阵
        /// </summary>
        double[,] m_pStreamNet;			//河网水系矩阵
        /// <summary>
        /// //河道Strahler等级
        /// </summary>
        public int[,] m_pStrahlerOrd;		//河道Strahler等级
        /// <summary>
        /// //河道最长长度计算
        /// </summary>
        double[,] m_pMaxFlowLength;		//河道最长长度计算
        /// <summary>
        /// 河道总长度计算
        /// </summary>
        double[,] m_pTotalFlowLength;	//河道总长度计算
        /// <summary>
        /// 流域边界矩阵
        /// </summary>
        double[,] m_pWatershed;			//流域边界矩阵
        /// <summary>
        /// 子流域栅格矩阵
        /// </summary>
        int[,] m_pSubWatershed;		//子流域栅格矩阵
        Cell FindedCell;
        private int PixelCount; //用于统计累积流量的变量
        public SortedDictionary<Cell, int> RiverNet = null;
        /// <summary>
        /// 包含所有与水体邻接的非水体像素的集合
        /// </summary>
        private SortedDictionary<Cell, int> TSet = new SortedDictionary<Cell, int>();
        //DEM栅格汇流参数
        public double m_pRouteCode;
        public   double[,] m_pRouteOrd;
        public double m_pRouteSequ;
        /// <summary>
        /// 汇流长度，即从某像素汇流到河口断的距离
        /// </summary>
        public double[,] m_pGridRouteLength;

        public double m_pGridRiverAvSlp;
        /// <summary>
        /// 地表径流汇流时间
        /// </summary>
        public double[,] m_pSurfQRouteTime;
        /// <summary>
        /// 壤中流汇流时间
        /// </summary>
        public double[,] m_pLatQRouteTime;
        /// <summary>
        /// 地下径流汇流时间
        /// </summary>
        public double[,] m_pBaseQRouteTime;

    }
    public struct Cell:IComparable<Cell>
    {
        public int row,col;//单元的行标x、列标y
        
        public Cell(int row, int col)
        {
            this.row = row;
            this.col = col;
        }
        public int CompareTo(Cell cell)
        {
            if (this.row < cell.row)
                return -1;
            else if (this.row > cell.row)
                return 1;
            else //==
            {
                if(this.col < cell.col)
                    return -1;
                if(this.col>cell.col)
                    return 1;
            }
            return 0;//二者完全相等
        }

    }
}
